%%************************************************************************************************************
% CPM(L=2,M=4,h=1/4) 的调制和维特比译码（对调制信号的差分信号译码）
% 经过相位倾斜后，一共有16个状态需要判断译码。
% 加入噪声比较分支度量的判决效果
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
close all;

%%% CPM中的由g(t)求积分q(t)
id=sqrt(-1);
pi=3.14159;
h=1/4;                                       % 调制指数  
M=4;                                         % 进制数  
% bit_number=30000;                            %540
bit_number=100;
symbol_number=bit_number/(log2(M));          % 码元数目
sample_number=8;                             % 每码元样点个数
sam=sample_number;
Tb=1/32000;                                  % 码元持续时间  ,rb=24k,tb=1/rb
fs=1/Tb*sam;                                 % 基带取样率

% data 产生0，1比特流
data=zeros(1,bit_number);
for i=1:bit_number
    temp=rand;
    if (temp<1/2)
        data(i)=0;
    elseif (temp<1) 
       data(i)=1;
   end
end;
%----------------------------------------------------------------------------------------------------------------------
% Ik={+1,-1,...+(M-1),-(M-1)} M=8  将比特流映射成M进制码元
Ik=zeros(1,symbol_number);
Uk=zeros(1,symbol_number);
for i=1:symbol_number
   if (data(2*i-1)==0&data(2*i)==0)
       Ik(i)=-3;
   elseif (data(2*i-1)==0&data(2*i)==1)
       Ik(i)=-1;
   elseif (data(2*i-1)==1&data(2*i)==1)
       Ik(i)=1;
   elseif (data(2*i-1)==1&data(2*i)==0)
       Ik(i)=3;
   end
end
    
% Ik=randi([0, 3], 1,symbol_number);
% for n=1:symbol_number
%     Ik(n)=Ik(n)*2-3;
% end
% Ik=[3	-1	3  -1	3	1	3	-3	-1	 1	 1	 3	 1	-3	-1 ...
%     -1	 1	1	1  -1	3	1	-1	 1	-1	-3	 1	-1	-1	-3 ...	
%     -1	-1 -3  -1  -1	1  -1	 1	 1	-3	 3	-3	 3	 1	 3 ...
%     -3	 1 -1	3  -1	3  -1	 1	 3	-1	-3	-3	 3	 3	 3 ...
%     -3	 1 -1  -1	3  -3  -1	-1	 1	-1	-1	-1	 3	-3	-1 ...
%      1	-1	1	-1	1	3  -3	-1	-3	 1	 3	 1	 1	-1	-1 ...
%     3	-3	3	-1	3	3	-1	-1	-1	3	3	-3	3	-1	3 ...	
%     -1	1	1	-3	3	3	3	3	-1	-1	3	-3	1	-3	-1 ...
%     3	-1	-3	1	3	-3	-1	-1	-1	-1	-3	-3	-1	-1	3 ...
%     -1	1	-3	-3	-1	1	1	1	-3	-3	-1	1	3	-1	1 ...
%     -1	1	-1	-1	1	-3	-1	3	-1	1	1	1	3	-3	3 ...	
%     -1	-3	1	3	-1	1	1	-3	1	-1	-3	-1	3	-1	1 ...
%     -3	3	-1	3	-3	3	3	1	3	3	-1	1	3	-1	-1 ...	
%     1	-1	-3	3	-1	1	3	1	-3	-1	-3	-1	1	1	3 ...
%     -1	1	-3	3	3	-3	-1	3	-3	1	-1	1	3	1	1 ...
%     3	3	3	-1	-3	3	1	3	3	1	-3	-3	-1	-3	3 ...
%     1	-3	-1	3	1	-3	-1	1	-1	3	1	-1	1	1	-3 ...	
%     1	1	-3	-1	1	-1	-1	-1	3	-1	1	-1	3	-3	-1];
%%%%%%%%------- M=4,L=2,h=0.25
L=2;
g=zeros(1,sam+1);
for i=0:L*sam
    t=i/fs;
    g(i+1)=1/(2*L*Tb)*(1-cos(2*pi*t/(L*Tb)));
end;
q=zeros(1,L*sam);
for i=1:L*sam
    for j=1:i
        q(i)=q(i)+g(j)/fs;
    end;
end;
qq=zeros(1,L*sam);
for i=1:L*sam
    qq(i)=2*pi*h*q(i);
end;
qtheta=1000*qq;
%figure(3);
%subplot(3,1,1);
%plot(g);
%subplot(3,1,2);
%plot(q);
%subplot(3,1,3);
%plot(qq);

seta=zeros(1,length(Ik)+1);
fii=zeros(1,length(Ik)*sam);
for j=2:length(Ik)+1
    seta(j+1)=seta(j)+pi*h*Ik(j-1);
end
for k=1:sam
    fii(k)=qq(k)*Ik(1);
    newfi(k)=fii(k)+pi*h*(M-1)*k/sam;
end
for j=2:length(Ik)  
    for i=1:sam
        fii((j-1)*sam+i)=seta(j)+qq(i)*Ik(j)+qq(sam+i)*Ik(j-1);
        newfi((j-1)*sam+i)=fii((j-1)*sam+i)+(pi*h*(M-1)*((j-1)*sam+i))/sam;%%%%倾斜相位
    end
end
for i=1:length(fii)
    ca(i)=cos(fii(i));                  % I_basewave
    sa(i)=sin(fii(i));                  % Q_basewave
end;
%  figure(1);
%  subplot(3,1,1);
%  plot(seta);
%  subplot(3,1,2);
%  plot(fii);
% subplot(3,1,3);
% plot(newfi);

% figure(2);
% subplot(2,1,1);
% plot(ca);
% subplot(2,1,2);
% plot(sa);

%%%%%%%%%%%%%%%%%%%%%%%%% 求发端基带信号差分 %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 验证发送是否正确 %%%%%%%%%%%%%%%%%%%%%%
% Sendrii=zeros(1,symbol_number*sample_number);
% Sendrqq=zeros(1,symbol_number*sample_number);
% for n=1:symbol_number
%  if n==1
%     for k=1:sam
  
%          Sendrii((n-1)*sam+k)=ca((n-1)*sam+k);
%          Sendrqq((n-1)*sam+k)=sa((n-1)*sam+k); %第一个码元的sam个采样点不作差分，直接送给差分值缓冲区

%     end 
%   else    
%     for k=1:sam
%         Sendrii((n-1)*sam+k)=ca((n-1)*sam+k)*ca((n-2)*sam+k)+sa((n-1)*sam+k)*sa((n-2)*sam+k);  %%cos(fk)=I(k)*I(k-1)+Q(k)*Q(k-1)
%         Sendrqq((n-1)*sam+k)=sa((n-1)*sam+k)*ca((n-2)*sam+k)-ca((n-1)*sam+k)*sa((n-2)*sam+k);  %%sin(fk)=Q(k)*I(k-1)-I(k)*Q(k-1)
%     end
%   end
% end
% figure(1);
% subplot(2,1,1);
% plot(Sendrii);
% subplot(2,1,2);
% plot(Sendrqq);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%*********************************基带加上噪声s+noise
ibase=awgn(ca,5,'measured');
qbase=awgn(sa,5,'measured');
%figure(3);
%subplot(2,1,1);
%plot(ibase);
%subplot(2,1,2);
%plot(qbase);
iwave=interp(ibase,4);               %%  96k/24k=4
qwave=interp(qbase,4);
%*****************%%%%%%%%%%%%%%%%%%%%%55 乘载波
% iwave=interp(ca,4);               %%  96k/24k=4
% qwave=interp(sa,4);

%  figure(4);
%  subplot(2,1,1);
%  plot(iwave);
%  subplot(2,1,2);
%  plot(qwave);

fc=128000;
Fs=fc*sam;             %% sample rate =fc*sam
t=0:1/Fs:(length(iwave)-1)*1/Fs;
isignal=iwave.*cos(2*pi*fc*t);
qsignal=qwave.*cos(2*pi*fc*t);
modsignal=iwave.*cos(2*pi*fc*t)-qwave.*sin(2*pi*fc*t);
mod=qwave.*sin(2*pi*fc*t)-iwave.*cos(2*pi*fc*t);

% figure(4);
% subplot(2,1,1);
% plot(isignal);
% subplot(2,1,2);
% plot(qsignal);

% figure(5);
% plot(mod);

%55555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555
%%%%%%%%%%%%%%%%%%% demodulate 
N=128;                                              % 滤波器的阶数为(N+1)  
F=[0,40000,45000,Fs/2]*2/Fs;
A=[1,1,0,0];
lpf=firls(N,F,A);
[amp_lpf,w]=freqz(lpf);
% plot(lpf);

% figure(1);
% grid on;
% plot(w/pi*(Fs/2),abs(amp_lpf),'b');                       % 低通滤波器频响特性 

%M=9;
fcc=fc-h*(M-1)/(2*Tb);                             %%%%倾斜相位解调
%fcc=85333;                                         %%%96k*8/9,倾斜相位一个周期采样9个点，对译码结果没有影响
t=0:1/Fs:(length(modsignal)-1)*1/Fs;
Ides=modsignal.*cos(2*pi*fcc*t)*2;
Ides=conv(Ides,lpf);
Ides=Ides(N/2+1:N/2+length(modsignal));   %%sa函数以N/2对称，卷积后只取数据与后面N/2个点的卷积，前面与之对称
Qdes=modsignal.*(-sin(2*pi*fcc*t)*2);
Qdes=conv(Qdes,lpf);
Qdes=Qdes(N/2+1:N/2+length(modsignal));

% figure(5);
% subplot(2,1,1);
% plot(Ides);
% subplot(2,1,2);
% plot(Qdes);

I_out=zeros(1,length(Ides)/4);            %% 抽取
Q_out=zeros(1,length(Ides)/4);
for i=1:length(I_out)
   I_out(i)=Ides(4*(i-1)+1);
   Q_out(i)=Qdes(4*(i-1)+1);
end;
% plot(ca)
% figure;
% hold on;
% plot(I_out,'r');


%ii=I_out(17:)
% figure(6);
% subplot(2,1,1);
% plot(I_out);
% subplot(2,1,2);
% plot(Q_out);

%－－－－－－－－－－－－****一个码元8个样点
% Icos=fix(I_out*0.95*32768);
% Qsin=fix(Q_out*0.95*32768);
% fid = fopen('d:\\Icos.txt','w');      %%%%倾斜相位后的基带信号－－I路
% for n =1:length(Icos)/sam
% fprintf(fid,' .word\t%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d\n',Icos(n*8-7),Icos(n*8-6),Icos(n*8-5),Icos(n*8-4),Icos(n*8-3),Icos(n*8-2),Icos(n*8-1),Icos(n*8));
% end
% fclose(fid);

% fid = fopen('d:\\Qsin.txt','w');     %%%%%倾斜相位后的基带信号－－Q路
% for n =1:length(Icos)/sam
% fprintf(fid,' .word\t%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d\n',Qsin(n*8-7),Qsin(n*8-6),Qsin(n*8-5),Qsin(n*8-4),Qsin(n*8-3),Qsin(n*8-2),Qsin(n*8-1),Qsin(n*8));
% end
% fclose(fid);

%%%%%%% I,Q倾斜后的基带数据存在一个缓冲区  %%%%%%%

% Icos=fix(I_out*0.95*32768);
% Isin=fix(Q_out*0.95*32768);
% fid=fopen('d:\\IQTitledData.txt','w');
% for n=1:length(Icos)/sam
%     fprintf(fid,' .word\t%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d\n',Icos(n*8-7),Isin(n*8-7),Icos(n*8-6),Isin(n*8-6),Icos(n*8-5),Isin(n*8-5),Icos(n*8-4),Isin(n*8-4),Icos(n*8-3),Isin(n*8-3),Icos(n*8-2),Isin(n*8-2),Icos(n*8-1),Isin(n*8-1),Icos(n*8),Isin(n*8));
% end
% fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

s_noise=I_out+Q_out*id;                       %% 基带信号--e^(j*seta)=cos(seta)+j*sin(seta)
r=s_noise;
len=length(r)/sam;
for n=1:len
 if n==1
    for k=1:sam
         rii((n-1)*sam+k)=I_out((n-1)*sam+k);
         rqq((n-1)*sam+k)=Q_out((n-1)*sam+k); %% 第一个码元的sam个采样点不作差分，直接送给差分值缓冲区
     end 
  else    
    for k=1:sam
        rii((n-1)*sam+k)=I_out((n-1)*sam+k)*I_out((n-2)*sam+k)+Q_out((n-1)*sam+k)*Q_out((n-2)*sam+k);  %% cos(fk)=I(k)*I(k-1)+Q(k)*Q(k-1)
        rqq((n-1)*sam+k)=Q_out((n-1)*sam+k)*I_out((n-2)*sam+k)-I_out((n-1)*sam+k)*Q_out((n-2)*sam+k);  %% sin(fk)=Q(k)*I(k-1)-I(k)*Q(k-1)
    end
  end
end

%  figure(10);
%  subplot(2,1,1);
%  plot(rii);
%  subplot(2,1,2);
%  plot(rqq);

%－－－－－－－－－－－－－－**一个码元16个样点的差分信号
%  Ide=fix(rii*0.95*32768);
%  Qde=fix(rqq*0.95*32768);

% fid = fopen('d:\\Ide.txt','w');      %%%%倾斜相位后的基带差分信号－－I路
% for n =1:length(Ide)/sam
% fprintf(fid,' .word\t%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d\n',Ide(n*8-7),Ide(n*8-6),Ide(n*8-5),Ide(n*8-4),Ide(n*8-3),Ide(n*8-2),Ide(n*8-1),Ide(n*8));
% end
% fclose(fid);

% fid = fopen('d:\\Qde.txt','w');     %%%%%倾斜相位后的基带差分信号－－Q路
% for n =1:length(Ide)/sam
% fprintf(fid,' .word\t%8d,%8d,%8d,%8d,%8d,%8d,%8d,%8d\n',Qde(n*8-7),Qde(n*8-6),Qde(n*8-5),Qde(n*8-4),Qde(n*8-3),Qde(n*8-2),Qde(n*8-1),Qde(n*8));
% end
% fclose(fid);

W=zeros(1,sam);
Ww=zeros(1,sam);
for i=1:sam
  W(i)=pi*h*(M-1)*i/sam-2*pi*h*(M-1)*(q(sam+i)+q(i))+(M-1)*pi*h;
   Ww(i)=pi*h*(M-1)*i/sam-(M-1)*(qq(sam+i)+qq(i))+(M-1)*pi*h;   %%qq(t)=q(t)*2*pi*h
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Uk 是倾斜相位的码元，(Uk Vn): Uk=(Ik+(M-1))/2; Vn=modp[U0+U1+U2+...+U(n-L)]
% Ik 是没有倾斜相位的码元，(Ik  THETAn)----即状态表示成(码元  相位)
% s_noise是倾斜相位后的调制信号，已经加过噪声
% r_T是差分的信号，就是s_noise(t)*s_noise(t-T)的共轭
% 状态数目N：h=m/p,m是偶数--N＝p*M^(L-1),M是奇数--N＝2*p*M^(L-1)；
% 奇偶时刻的状态不同。但当相位倾斜后，状态不分奇偶时刻，所以h=1/4,M=4,L=2时，N＝p*M^(L-1)＝16
% 转移状态矩阵transtate[4,4]。第i行表示第i个状态是从前一时刻的哪个状态转移过来的。(Un_1 Vn)表示一个状态，Vn=Vn_1+Un_2
%****64个状态转移时有模8现象：即第1,9,17,25,33,41,49,57状态转移和第1个状态的转移是一样的。
%****第8,16,24,32,40,48,56,64状态转移和第8个状态的转移是一样的，在除以8取余＝0时要加上8。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% viterbi 检测
%----------------------------16 states
%state_number=64;                             %% 状态数目：
%delay=10;                                    %% 维特比译码的码元段长
%len=length(r)/sam;  %% sam=5(一个码元的样点数目),len 是码元个数
%transtate=[1 16 23 30 37 44 51 58;
           % 2 9 24 31 38 45 52 59;
           % 3 10 17 32 39 46 53 60;
           % 4 11 18 25 40 47 54 61;
           % 5 12 19 26 33 48 55 62;
           % 6 13 20 27 34 41 56 63;
           % 7 14 21 28 35 42 49 64;
           % 8 15 22 29 36 43 50 57];  
%metric_state=zeros(state_number,delay);     %% 矩阵(64,10)，存储的是最大度量和的值，列表示时刻
%survivor_state=zeros(state_number,delay);   %%幸存状态(64,10)，列表示时刻，行表示状态
%decision=zeros(1,len);                      %% 判断出来的码元－－Uk
%Ik_decision=zeros(1,symbol_number);
%testcos=zeros(1,512*5);
%testsin=zeros(1,512*5);

state_number=4;
delay=10;
len=length(s_noise)/sample_number;
transtate=[1,2,3,4]; 
estimated_state_old=[1,5,9,13];  %%除了1，其他可以是简化前的任意状态，只要是在简化范围内就行。
estimated_state_old2=[1,5,9,13];  %%就是上二时刻的状态，是简化前的4个
estimated_state_new=zeros(1,state_number);     %%当前的状态，是简化前的4个。
metric_state=ones(state_number,delay);         %%元素取值范围是1--4，矩阵是4*10
survivor_state=zeros(state_number,delay);      %%元素取值范围为1--4，矩阵是4*10
decision=zeros(1,len);
Ik_decision=zeros(1,symbol_number);

% qwu=zeros(16,len);
%--------------------------------------------------------------------------
for n=1:len     %% n是码元个数
   %*******************************************************************************************************************
  if n==1
    %-------------------------------------------
    Un_2=0;
    Vn_1=0;
    for i=1:state_number          %%  i是状态数目4
        Un=i-1;                   %%  i是当前状态
        for j=1:4
            last_state=estimated_state_old(transtate(j)); %% 还是简化前的状态编号，
            Un_1=fix((last_state-1)/4);
%             ki=4*(i-1)+j;
%             qwu(ki,n)=Un_1;
            Vn=rem((last_state-1),4);
            for k=1:sam   
                 reffi1(k)=2*pi*h*Vn+4*pi*h*Un*q(k)+4*pi*h*Un_1*q(sam+k);             %% 当前时刻的角度
                 reffi2(k)=2*pi*h*Vn_1+4*pi*h*Un_1*q(k)+4*pi*h*Un_2*q(sam+k)+W(k);    %% 前一时刻的角度
                 I_ref(k)=cos(reffi1(k))*cos(reffi2(k))+sin(reffi1(k))*sin(reffi2(k));%% 参考相位的差分余弦值
                 Q_ref(k)=sin(reffi1(k))*cos(reffi2(k))-cos(reffi1(k))*sin(reffi2(k));%% 参考相位的差分正弦值
                      mi(k)=rii((n-1)*sam+k)*I_ref(k);
                      mq(k)=rqq((n-1)*sam+k)*Q_ref(k);
            end
            metric=sum(mi+mq);
            if metric>metric_state(i,n)
                metric_state(i,n)=metric;
                survivor_state(i,n)=transtate(j);  
                Vn=rem((Un_1+Vn),4);               %% Vn就是V(n+1)=rem((Un_1+Vn),4),刷新
                Un_1=Un;            
                estimated_state_new(i)=Un_1*4+Vn+1;%% 此时刻的状态在下一时刻就是last_state,所以放在estimated_state_new中间
           end
       end  % end j=4
    end     % end i=4
   %*******************************************************************************************************************    
  elseif n<=delay
    %---------------------------------------------------------------------
    for i=1:state_number
        Un=i-1;      %%%% i 是当前状态
        for j=1:4
            last_state=estimated_state_old(transtate(j));
            last2_state=estimated_state_old2(survivor_state(transtate(j),n-1));
            Un_1=fix((last_state-1)/4);
%             ki=4*(i-1)+j;
%             qwu(ki,n)=Un_1;
            Vn=rem((last_state-1),4);
            if last2_state==0
                Un_2=0;
                Vn_1=0;
            else Un_2=fix((last2_state-1)/4);
                 Vn_1=rem((last2_state-1),4);
            end
            for k=1:sam
                 reffi1(k)=2*pi*h*Vn+4*pi*h*Un*q(k)+4*pi*h*Un_1*q(sam+k);     %% 当前时刻的角度
                 reffi2(k)=2*pi*h*Vn_1+4*pi*h*Un_1*q(k)+4*pi*h*Un_2*q(sam+k); %% 前一时刻的角度
                 I_ref(k)=cos(reffi1(k))*cos(reffi2(k))+sin(reffi1(k))*sin(reffi2(k));  %% 参考相位的差分余弦值
                 Q_ref(k)=sin(reffi1(k))*cos(reffi2(k))-cos(reffi1(k))*sin(reffi2(k));  %% 参考相位的差分正弦值
                        mi(k)=rii((n-1)*sam+k)*I_ref(k);
                        mq(k)=rqq((n-1)*sam+k)*Q_ref(k);
            end 
            metric=sum(mi+mq);
           if metric+metric_state(transtate(j),n-1)>metric_state(i,n)
               metric_state(i,n)=metric+metric_state(transtate(j),n-1);   %% 保存最大的度量和
               survivor_state(i,n)=transtate(j);     %% transtate(current_state_number_moduloM,j);  %% 保存最大可能的上一时刻的状态编号
               Vn=rem((Un_1+Vn),4); 
               Un_1=Un;
               estimated_state_new(i)=Un_1*4+Vn+1;
           end
       end  % end j=4
    end     % end i=4
  %******************************************************************************************************************* 
  else        %%%% n>delay 的情况，delay相当于一个长度10的窗，在移动。大于10后，要刷新delay
      
    for i=1:state_number
        for t=1:delay-1     %%每一列顺序前移，第十列等待刷新，被计算出来的新值代替
           metric_state(i,t)=metric_state(i,t+1);
           survivor_state(i,t)=survivor_state(i,t+1);
        end
    end
    %metric_state(:,delay)=zeros(state_number,1);
    %survivor_state(:,delay)=zeros(state_number,1);
    %------------------------------------------
      for i=1:state_number
           Un=i-1;
          for j=1:4
             last_state=estimated_state_old(transtate(j));
              last2_state=estimated_state_old2(survivor_state(transtate(j),delay-1));
              Un_1=fix((last_state-1)/4);
%               ki=4*(i-1)+j;
%             qwu(ki,n)=Un_1;
              Vn=rem((last_state-1),4);
              Un_2=fix((last2_state-1)/4);
              Vn_1=rem((last2_state-1),4);
              for k=1:sam
                 reffi1(k)=2*pi*h*Vn+4*pi*h*Un*q(k)+4*pi*h*Un_1*q(sam+k);     %%当前时刻的角度
                 reffi2(k)=2*pi*h*Vn_1+4*pi*h*Un_1*q(k)+4*pi*h*Un_2*q(sam+k); %%前一时刻的角度
                 I_ref(k)=cos(reffi1(k))*cos(reffi2(k))+sin(reffi1(k))*sin(reffi2(k));  %%参考相位的差分余弦值
                 Q_ref(k)=sin(reffi1(k))*cos(reffi2(k))-cos(reffi1(k))*sin(reffi2(k));  %%参考相位的差分正弦值
                 mi(k)=rii((n-1)*sam+k)*I_ref(k);
                 mq(k)=rqq((n-1)*sam+k)*Q_ref(k);
             end 
             metric=sum(mi+mq);
              if  metric+metric_state(transtate(1,j),delay-1)>metric_state(i,delay)
                 
                  metric_state(i,delay)=metric+metric_state(transtate(j),delay-1);
                  survivor_state(i,delay)=transtate(j);   %transtate(current_state_number_moduloM,j);
                  Vn=rem((Un_1+Vn),4);
                  Un_1=Un;
                  estimated_state_new(i)=Un_1*4+Vn+1;   
              end
          end  % end j=4 
    end        % end i=4
end            % else end 
  %%******************************************************************************************************************
  %% 每到一个delay的长度，就根据度量最大的值回溯判断出码元 
    %--------------------------------------------------------------------
    if (n>=delay)&(n<len)     %% 从n=delay开始回溯判断
        max=1;
        for i=2:state_number
           if metric_state(i,delay)>metric_state(max,delay)
              max=i;
           end
        end
        survivor=max;      %% 在delay列中最大度量值对应的行号.survivor_state中第i行就是第i个状态在不同时刻的前一状态，列表示时刻
        for t=delay:-1:2   %% 回溯到第2列－－第2列就是第一时刻的状态
            survivor=survivor_state(survivor,t);              %% 第j列的元素就是在第(j-1)时刻状态
        end
        decision(n-delay+1)=survivor-1;                       %% Uk就是倾斜后的码元0
        Ik_decision(n-delay+1)=2*decision(n-delay+1)-(M-1);   %% Uk=(Ik+(M-1))/2
       
    %--------------------------------------------------------------------
    elseif n==len     %% 最后10个码元就由metric矩阵判断
        max=1;
        for i=2:state_number
           if metric_state(i,delay)>metric_state(max,delay)
              max=i;
           end
        end
        survivor=max;  %% 在delay列中最大度量值对应的行号.
        for t=delay:-1:1
           decision(n-delay+t)=survivor-1;
           Ik_decision(n-delay+t)=2*decision(n-delay+t)-(M-1);     %% Uk=(Ik+(M-1))/2
           survivor=survivor_state(survivor,t);
        end
   end % for elseif
   %--------------------------------------------------------------------
   estimated_state_old2=estimated_state_old;
   estimated_state_old=estimated_state_new;  
end    % for n=1:len end
newmetr=metric_state';
for i=1:10
    for j=1:4
    duliang((i-1)*8+j)=newmetr(i,j);
    end
end
%figure(19);
%plot(duliang);
%figure(20);
%subplot(2,1,1);
%plot(Uk);
%plot(Ik);
%subplot(2,1,2);
%plot(Ik_decision);
figure;
plot(Ik_decision-Ik);
for i=1:symbol_number 
    erro(i)=Ik_decision(i)-Ik(i);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:symbol_number
    if Ik_decision(i)==-3
        data_decision(2*i-1)=0;
        data_decision(2*i)=0;
    elseif Ik_decision(i)==-1
        data_decision(2*i-1)=0;
        data_decision(2*i)=1;
    elseif Ik_decision(i)==1
        data_decision(2*i-1)=1;
        data_decision(2*i)=1;
    elseif Ik_decision(i)==3
        data_decision(2*i-1)=1;
        data_decision(2*i)=0;
    end
end
     
[symbolerror_num,symbolerror_rat]=symErr(Ik,Ik_decision)      % symbol error rate
% [biterror_num,biterror_rat]=symErr(data,data_decision)           % bit error rate
metric_state=zeros(state_number,delay);         %%元素取值范围是1--4，矩阵是4*10
survivor_state=zeros(state_number,delay);      %%元素取值范围为1--4，矩阵是4*10
decision=zeros(1,len);
Ik_decision=zeros(1,symbol_number);
tab=zeros(128,8);
for i_0=0:3
    for i_2=0:3
        for i_1=0:3
            li=i_0*16+i_2*4+i_1;
            for k=1:4
            tab(2*li+1,2*k-1)=cos(2*pi*h*i_2+4*pi*h*(i_0-i_1)*q(k)+4*pi*h*q(sam+k)*(i_1-i_2));
            tab(2*li+1,2*k)=sin(2*pi*h*i_2+4*pi*h*(i_0-i_1)*q(k)+4*pi*h*q(sam+k)*(i_1-i_2));
            k1=k+4;
            tab(2*li+2,2*k-1)=cos(2*pi*h*i_2+4*pi*h*(i_0-i_1)*q(k1)+4*pi*h*q(sam+k1)*(i_1-i_2));
            tab(2*li+2,2*k)=sin(2*pi*h*i_2+4*pi*h*(i_0-i_1)*q(k1)+4*pi*h*q(sam+k1)*(i_1-i_2));
        end
        end
    end
end
%--------------------------------------------------------------------------
 icf=1;%
    Vn_1=0;Un_2=zeros(1,4);
    for i=1:state_number          %%  i是状态数目4
        Un=i-1;                   %%  i是当前状态
        for j=1:4
%            last_state=estimated_state_old(transtate(j)+1); %% 还是简化前的状态编号，
            Un_1=j-1;%fix((last_state-1)/4);
            
           % Vn=rem((last_state-1),4);
             li=16*Un+4*Un_2(j)+Un_1;
                      for k=1:4
                   I_ref(k)=tab(2*li+1,2*k-1);
                   Q_ref(k)=tab(2*li+1,2*k);
                   I_ref(k+4)=tab(2*li+2,2*k-1);
                   Q_ref(k+4)=tab(2*li+2,2*k);
                      end
             for k=1:icf:sam   
%                  reffi1(k)=4*pi*h*Un*q(k)+4*pi*h*Un_1*q(sam+k);             %% 当前时刻的角度
%                  reffi2(k)=4*pi*h*Un_1*q(k)+W(k);    %% 前一时刻的角度
%                  I_ref(k)=cos(reffi1(k))*cos(reffi2(k))+sin(reffi1(k))*sin(reffi2(k));%% 参考相位的差分余弦值
%                  Q_ref(k)=sin(reffi1(k))*cos(reffi2(k))-cos(reffi1(k))*sin(reffi2(k));%% 参考相位的差分正弦值
                      mi(k)=rii((n-1)*sam+k)*I_ref(k);
                      mq(k)=rqq((n-1)*sam+k)*Q_ref(k);
            end
                metric=sum(mi+mq);
             if metric>metric_state(i,1)%
                metric_state(i,1)=metric;
                survivor_state(i,1)=j; %transtate(j)
%                Vn=rem((Un_1+Vn),4);               %% Vn就是V(n+1)=rem((Un_1+Vn),4),刷新
                Un_22(i)=j-1;
              %  Un_1=Un;            
              %  estimated_state_new(i)=Un_1*4+Vn+1;%% 此时刻的状态在下一时刻就是last_state,所以放在estimated_state_new中间
            end%
       end  % end j=4
    end     % end i=4
    % estimated_state_old2=estimated_state_old;
    % estimated_state_old=estimated_state_new;  
   %*******************************************************************************************************************    
   Un_2=Un_22;
   %elseif n<=delay
   I_ref=zeros(1,8);
    Q_ref=zeros(1,8);
   for n=2:delay-1
    %---------------------------------------------------------------------
      for i=1:state_number
           Un=i-1;      %%%% i 是当前状态
           for j=1:4
               Un_1=j-1;
                  li=16*Un+4*Un_2(j)+Un_1;
                      for k=1:4
                   I_ref(k)=tab(2*li+1,2*k-1);
                   Q_ref(k)=tab(2*li+1,2*k);
                   I_ref(k+4)=tab(2*li+2,2*k-1);
                   Q_ref(k+4)=tab(2*li+2,2*k);
                      end
%                    I_ref=tab(2*li+1,1:8);
%                    Q_ref=tab(2*li+2,1:8);
               for k=1:icf:sam
                   %reffi1(k)=2*pi*h*Un_2(j)+4*pi*h*(Un-Un_1)*q(k)+4*pi*h*(Un_1-Un_2(j))*q(sam+k);
                   %I_ref(k)=cos(reffi1(k));%reffi1(k)-reffi2(k)
                     % Q_ref(k)=sin(reffi1(k));%reffi1(k)-reffi2(k)
                   mi(k)=rii((n-1)*sam+k)*I_ref(k);
                   mq(k)=rqq((n-1)*sam+k)*Q_ref(k);
               end 
               metric=sum(mi+mq);
              if metric+metric_state(j,n-1)>metric_state(i,n)
                   metric_state(i,n)=metric+metric_state(j,n-1);   %% 保存最大的度量和transtate(j)+1
                   survivor_state(i,n)=j;     %% transtate(current_state_number_moduloM,j);  %% 保存最大可能的上一时刻的状态编号transtate(j)
%                   Vn=rem((Un_1+Vn),4); 
                   Un_22(i)=j-1;
                   %j-1
                 %  Un_1=Un;
%                   estimated_state_new(i)=Un_1*4+Vn+1;
              end
           end   % end j=4
           
    end     % end i=4
     %estimated_state_old2=estimated_state_old;
   %  estimated_state_old=estimated_state_new;  
     Un_2=Un_22;
 end
  %******************************************************************************************************************* 
  %else        %%%% n>delay 的情况，delay相当于一个长度10的窗，在移动。大于10后，要刷新delay
  for n=delay:len   
         for i=1:state_number
              Un=i-1;
              for j=1:4
                   Un_1=j-1;
                    li=16*Un+4*Un_2(j)+Un_1;
                    for k=1:4
                   I_ref(k)=tab(2*li+1,2*k-1);
                   Q_ref(k)=tab(2*li+1,2*k);
                   I_ref(k+4)=tab(2*li+2,2*k-1);
                   Q_ref(k+4)=tab(2*li+2,2*k);
               end
                   for k=1:icf:sam
%                       reffi1(k)=2*pi*h*Un_2(j)+4*pi*h*(Un-Un_1)*q(k)+4*pi*h*(Un_1-Un_2(j))*q(sam+k);%
%                       I_ref(k)=cos(reffi1(k));%reffi1(k)-reffi2(k));%
%                       Q_ref(k)=sin(reffi1(k));%reffi1(k)-reffi2(k));% 
                     mi(k)=rii((n-1)*sam+k)*I_ref(k);
                      mq(k)=rqq((n-1)*sam+k)*Q_ref(k);
                   end 
                   metric=sum(mi+mq);
                   if  metric+metric_state(j,delay-1)>metric_state(i,delay)
                       metric_state(i,delay)=metric+metric_state(j,delay-1);
                       survivor_state(i,delay)=j;   %transtate(current_state_number_moduloM,j);transtate(j)
                       Un_22(i)=j-1;
                  end
              end  % end j=4 
        end        % end i=4
  
  %%******************************************************************************************************************
  %% 每到一个delay的长度，就根据度量最大的值回溯判断出码元 
    %--------------------------------------------------------------------
           max=1;
         for i=2:state_number
             if metric_state(i,delay)>metric_state(max,delay)
                  max=i;
             end
        end
        survivor=max;      %% 在delay列中最大度量值对应的行号.survivor_state中第i行就是第i个状态在不同时刻的前一状态，列表示时刻
        for t=delay:-1:2   %% 回溯到第2列－－第2列就是第一时刻的状态
            survivor=survivor_state(survivor,t); 
            survivor=survivor;%% 第j列的元素就是在第(j-1)时刻状态+1
        end
        decision(n-delay+1)=survivor-1;                       %% Uk就是倾斜后的码元0
        Ik_decision(n-delay+1)=2*decision(n-delay+1)-(M-1);   %% Uk=(Ik+(M-1))/2
         Un_2=Un_22;
        estimated_state_old2=estimated_state_old;
        estimated_state_old=estimated_state_new;  
        for i=1:state_number
              for t=1:delay-1     %%每一列顺序前移，第十列等待刷新，被计算出来的新值代替
                   metric_state(i,t)=metric_state(i,t+1);
                   survivor_state(i,t)=survivor_state(i,t+1);
              end
         end
end   
    %--------------------------------------------------------------------
        max=1;
        for i=2:state_number
           if metric_state(i,delay)>metric_state(max,delay)
              max=i;
           end
        end
        Ik_decision(n)=2*(max-1)-(M-1);
        survivor=max;  %% 在delay列中最大度量值对应的行号.
        for t=delay-1:-1:1
            survivor=survivor_state(survivor,t);
           decision(n-delay+t)=survivor-1;
           Ik_decision(n-delay+t)=2*decision(n-delay+t)-(M-1);     %% Uk=(Ik+(M-1))/2
          % survivor=survivor+1;
        end
        
% end % for elseif
%    %--------------------------------------------------------------------
%    estimated_state_old2=estimated_state_old;
%    estimated_state_old=estimated_state_new;  
% end    % for n=1:len end
% newmetr=metric_state';
% for i=1:10
%     for j=1:4
%     duliang((i-1)*8+j)=newmetr(i,j);
%     end
% end
%figure(19);
%plot(duliang);
%figure(20);
%subplot(2,1,1);
%plot(Uk);
%plot(Ik);
%subplot(2,1,2);
figure;
plot(Ik_decision-Ik);
for i=1:symbol_number 
    erro(i)=Ik_decision(i)-Ik(i);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:symbol_number
    if Ik_decision(i)==-3
        data_decision(2*i-1)=0;
        data_decision(2*i)=0;
    elseif Ik_decision(i)==-1
        data_decision(2*i-1)=0;
        data_decision(2*i)=1;
    elseif Ik_decision(i)==1
        data_decision(2*i-1)=1;
        data_decision(2*i)=1;
    elseif Ik_decision(i)==3
        data_decision(2*i-1)=1;
        data_decision(2*i)=0;
    end
end
     
[symbolerror_num,symbolerror_rat]=symErr(Ik,Ik_decision)      % symbol error rate